//888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//	File		 :	AnnotateSequence.cpp
//	Contents	 :	Source file for functions annotating protein sequences
//						links functionality from parseMSF.h and FormattedOutput.h
//					Assumes only A-Z upper case letters represent residues
//					Uses Eiffel-type pre/postcondition assertions for more fragile code
//	Creation Date:	23/12/1999
//	Last Modified:	24/12/1999	lg
//	tabs width	 :	every 4 characters
//
//888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#include <stdexcept>
#include <vector>
#include <string>
#include <algorithm>
#include <cassert>

#include "annotatesequence.h"
#include "formattedoutput.h"
#include "sequence.h"
#include "parsefasta.h"
#include "parsemsf.h"
#include "parseclustal.h"
#include "parsephd.h"

using std::string;
using std::vector;

#define GAP_SYMBOL (unsigned char)26
#define XXX_SYMBOL (unsigned char)27

//_______________________________________________________________________________________
//
//	ResidueTypeArraysToSets
//
//		Description:
//			Residue Groups are strings each of whose characters represents an amino acid
//			 type.
//			For performance reasons, this function converts these strings into sets.
//
//		Implementation:
//			The sets are thirty-two bytes in length (>26 alphabetic characters)
//			The presence of each residue is flagged by turnning on all the bits
//			 in the corresponding byte.
//
//		Last modified: lg 24/12/1999
//
//_______________________________________________________________________________________
namespace
{
    void ResidueTypeArraysToSets(const CStrVec& arrResidueType, CSetsVec& residueGroups)
    {
        residueGroups.clear();
        residueGroups.reserve(arrResidueType.size());
        for (unsigned i = 0; i < arrResidueType.size(); ++i)
        {
            residueGroups.push_back(CSet32Short());
            residueGroups[i].Clear();
            for (int j = arrResidueType[i].size() - 1; j >= 0; --j)
            {
                // characters represent aa and should only be capital letters
                assert(isupper(arrResidueType[i][j]));
                residueGroups[i].SetFlag(arrResidueType[i][j] - 'A');
            }
        }

        // POSTCONDITION
        assert(arrResidueType.size() == residueGroups.size());
        // POSTCONDITIONS END
    }

    void CreateLookupTable(const std::string& invalidAAs, char* lookup, bool lowercase, bool ignoreGaps)
    {
        memset(lookup, GAP_SYMBOL, 256);
        for (char i = 'A'; i <= 'Z'; ++i)
            lookup[i] = i - 'A';
        for (int i = invalidAAs.size() - 1; i >= 0; --i)
            lookup[invalidAAs[i]] = GAP_SYMBOL;
        if (lowercase)
        {
            for (int i = invalidAAs.size() - 1; i >= 0; --i)
                lookup[invalidAAs[i] + ('a'-'A')] = GAP_SYMBOL;
            for (char i = 'a'; i <= 'z'; ++i)
                lookup[i] = i - 'a';
        }
		if (ignoreGaps)
		{
			if (lowercase)
				lookup['x'] = XXX_SYMBOL;
			lookup['X'] = XXX_SYMBOL;
		}
    }
}


//_______________________________________________________________________________________
//
//	CreateAnnotation
//
//		Description:
//			Takes a list of consensus types, a vector of residue sequences and
//				a consensus threshold.
//			Returns data in a CSequence format indicating which consensus residues
//				should be formatted in a particular way
///
//		Inputs:
//			residueGroups		: Sets of residue types (see above)
//			sequence			: Aligned amino acid sequences from different proteins
//			consensusPermille	: The threshold determining the presence of consensus
//									residues in tenths of a percent
//
//		Outputs:
//			sequence			: Aligned amino acid sequences from different proteins
//			formatSequence		: list of formats corresponding to the residues at the
//									various positions in sequence
//
//
//		Implementation:
//			Extensive use is made of set operations. Though the efficiency gains may not
//			 be very significant here, the logic can be reused where larger cross
//			 matches results in a hyperbolic explosion in the number of compares.
//			Naive code: compares = protein counts x residue group counts
//			Set code  : compares = residue group counts ( + overhead of set operations)
//
//		Last modified: lg 24/12/1999
//
//_______________________________________________________________________________________
void CreateAnnotation(	//inputs
    const CStrVec& aminoAcidSets,
    const std::string& markers,
    const std::string& invalidAAs,
    const CAnnotateOptions& opt,
	bool  chkAAIdentity,
	bool  allowLowerCase,
    // outputs
    CSequence& sequence,
    CSequence& formatSequence)
{
    //88888888888888888888888888888888888888888888888888888888888888888888888888888888888
	// PRECONDITIONS
	assert(sequence.IS_VALID());
	// Because of the use of byte sized sets,
	//	this function can only handle a max of 65536 proteins
	assert(sequence.Count() < 65536);
	// can there be a consensus < 50. This would allow > 1 consensus-on-1-residue?
	assert(opt.m_ConsensusThreshold >= 500);
    // PRECONDITIONS END
    //88888888888888888888888888888888888888888888888888888888888888888888888888888888888

	CSetsVec residueGroups;
	ResidueTypeArraysToSets(aminoAcidSets, residueGroups);

	// create lookup table for speed
	char lookup[256];
	CreateLookupTable(invalidAAs, lookup, allowLowerCase, opt.m_XForUnknown);


	// allocate room for the formats
	formatSequence.Len = sequence.Len;
	for (unsigned i = 0; i < sequence.Count() ; ++i)
		formatSequence.Data.push_back(std::string(sequence.Len, '\0'));
	vector<string>(sequence.Count()).swap(formatSequence.Names);
	
	// allocate room for the consensus line
	sequence.ExtraNames.insert(sequence.ExtraNames.begin(), "Markers");
	sequence.ExtraData.insert(sequence.ExtraData.begin(), std::string(sequence.Len,
																			markers[0]));


	CSet32Short aaTable, resultSet;

	// array of all the residues at this position
	vector<char> alignedAA(sequence.Count(), ' ');

	// use last format for single residue matches
	char singleFormat = (char)(residueGroups.size() + 1);

	// Move down amino acid sequence one position at a time
	for (unsigned pos = 0; pos < sequence.Len; ++pos)
	{
		//
		// Tabulate a set of all the different residues at this position
		//
		aaTable.Clear();
		for (int j = sequence.Count() - 1; j >= 0; --j)
			aaTable.IncFlag(alignedAA[j] = lookup[(unsigned char)sequence.Data[j][pos]]);

		//
		// Calculate consensus threshold
		//
		// N.B. The sequence consensus threshold is the next whole number of sequences
		//			greater than the real number calculated
		// i.e. not round() but ceiling()
		unsigned consensusThreshold = sequence.Count();
		if (!opt.m_IncludeBlanks)
			consensusThreshold -= aaTable.IsFlagged(GAP_SYMBOL);
		if (opt.m_XForUnknown)
			consensusThreshold -= aaTable.IsFlagged(XXX_SYMBOL);
		consensusThreshold = (consensusThreshold * opt.m_ConsensusThreshold + 999) / 1000;
		//
		// No consensus possible: all positions either gaps or XForUnknowns
		//
		if (!consensusThreshold)
			continue;


		//
		// consensus on just one residue (identity)
		//
		if (chkAAIdentity							&&
			aaTable.GetMaxIndex() < GAP_SYMBOL		&&
			aaTable.GetMaxCnt() >= consensusThreshold)
		{
			char maxResidue = (char)aaTable.GetMaxIndex();
			
			// marker is the consensus residue
			sequence.ExtraData[0][pos] = (char)('A' + maxResidue);

			//	flag all matching residues at this position for highlighting
			for (int i = sequence.Count() - 1; i >= 0; --i)
				if(maxResidue == alignedAA[i])
					formatSequence.Data[i][pos] = singleFormat;
			continue;
		}

		
		
		
		// Check residue groupings for consensus at this position
		//	i.e. if table of residues (aaTable) matches are over the consensus threshold
		// The set AND operator removes all residues not in the particular residue group
		unsigned consensus;
		for (consensus = 0; consensus < residueGroups.size(); ++consensus)
		{
			resultSet.AND(residueGroups[consensus], aaTable);
			if (resultSet.Total() >= consensusThreshold)
				break;

		}

		// consensus is the index of the aminoacid groups.

		if (consensus < residueGroups.size())
		{
			
			// Add 1 to skip over the default (null) format
			char format = (char)(consensus+1);
			
			// add marker
			sequence.ExtraData[0][pos] = markers[format];

			//	flag all consensus residues at this position for highlighting
			for (unsigned i = 0; i < sequence.Count(); ++i)
				if(residueGroups[consensus].IsFlagged(alignedAA[i]))
					formatSequence.Data[i][pos] = format;
		}
	}

    //88888888888888888888888888888888888888888888888888888888888888888888888888888888888
    // POSTCONDITIONS
	assert(formatSequence.IS_VALID());
	assert(sequence.IS_VALID());
    // POSTCONDITIONS END
    //88888888888888888888888888888888888888888888888888888888888888888888888888888888888
}


//_______________________________________________________________________________________
//
//	OutputAnnotatedSequence
//
//		Description:
//			Takes the lists of amino acids and drives the OutputEngine to
//				produce annotated sequences.
//			Change the output type by passing in a different OutputEngine
//			Two overloaded versions are provided, the 2nd of which is optimized
//				for when the number of protein sequences < 32.
//
//		Inputs:
//			sequence			: Aligned amino acid sequences from different proteins
//			formatSequence		: To hold the format (indices) of corresponding sequences
//									at corresponding positions.
//			wrap				: Long alignments wrap around every "wrap" residues
//									To avoid wrapping, set wrap = sequence.Len
//			columnW				: Divide annotation into columns of columnW width
//									To avoid columns, set columnW = sequence.Len
//			OutputEngine		: Use this object to produce annotated output in
//									different file formats
//
//		Last modified: lg 11/04/2000
//
//_______________________________________________________________________________________
void OutputAnnotatedSequence(	const CSequence& sequence,
							    const CSequence& formatSequence,
							    unsigned wrap,
							    unsigned columnW,
							    CFormattedOutput& OutputEngine)
{
    // PRECONDITIONS
	assert(sequence.IS_VALID());
	assert(formatSequence.IS_VALID());
	assert(wrap > 0 && wrap <= sequence.Len);
	assert(columnW > 0 && wrap <= sequence.Len);
    // PRECONDITIONS END


	for (unsigned beg = 0; beg < sequence.Len; beg += wrap)
	{
		unsigned end = std::min(sequence.Len, beg + wrap);
		for (unsigned i = 0; i < sequence.Count(); ++i)
		{
			//
			// start new line
			//
			if (i)
				OutputEngine.AddNewLine();
			
			//
			// add sequence name in default format followed by a single space
			//
			int currentFormat = 0;
			OutputEngine.SetFormat(0);
			OutputEngine.AddText(sequence.Names[i]);
			OutputEngine.AddText(' ');
			
			
			//
			// add sequence residues with optional column breaks
			//
			for (unsigned j = beg; j < end; ++j)
			{
				// columns breaks
				if (!((j - beg) % columnW))
				{
					if (0 != currentFormat)
						OutputEngine.SetFormat(currentFormat = 0);
					OutputEngine.AddText(' ');
				}
				
				// residues
				if (formatSequence.Data[i][j] != currentFormat)
                    OutputEngine.SetFormat(currentFormat = formatSequence.Data[i][j]);
				OutputEngine.AddText(sequence.Data[i][j]);
			}
		}
		
		//
		// Add extra line to separate consensus and other extra data lines
		//
		OutputEngine.AddNewLineAndHalf();
		OutputEngine.SetFormat(0);

		//
		// Add extra lines, e.g. for consensus, 2ndary structure
		//
		for (unsigned i = 0; i < sequence.ExtraCount(); ++i)
		{
			//
			// add extra names in default format followed by a single space
			//
			OutputEngine.AddText(sequence.ExtraNames[i]);
			OutputEngine.AddText(' ');
			
			//
			// Extra data with optional column breaks
			//
			for (unsigned j = beg; j < end; ++j)
			{
				// columns breaks
				if (!((j - beg) % columnW))
					OutputEngine.AddText(' ');
				
				// extra data
				OutputEngine.AddText(sequence.ExtraData[i][j]);
			}
            //			OutputEngine.AddText(sequence.ExtraData[i], beg, end - beg);
			OutputEngine.AddNewLine();
		}

		// add one line between wraps
		OutputEngine.AddNewLine();
	}


}

//_______________________________________________________________________________________
//
//	GuessAlignmentFileFormat
//
//		Guess the data format in an alignment buffer based on the first 2000 characters
//
//_______________________________________________________________________________________
CSequence::eFormatType GuessAlignmentFileFormat(char* buf, unsigned bufLen)
{
	assert(buf && bufLen);

	// skip past initial whitespace
	char* searchStr = buf;
	while (*searchStr && isspace(*searchStr))
		++searchStr;

//	CSequence::eFormatType formatType = CSequence::eUnknown;
	CSequence::eFormatType formatType = CSequence::eMSF;
	if (*searchStr == '>')
		formatType = CSequence::eFASTA;
	else if (searchStr + 7 < buf + bufLen && strncmp("CLUSTAL", searchStr, 7) == 0)
		formatType = CSequence::eCLUSTAL;
	else if (searchStr + 3 < buf + bufLen && strncmp("DSC", searchStr, 3) == 0)
		formatType = CSequence::eDSC;
	else
	{
		// limit search somehow
		char * endSearch = 0;
		char originalChar = 0;
		if (searchStr + 2000 < buf + bufLen)
		{
			endSearch = searchStr + 2000;
			originalChar = *endSearch;
			*endSearch = 0;
		}

		//search for token identifiers for the two file types
		
		if (char* phdstr = strstr(searchStr, "PHDsec"))
		{
			if (strstr(phdstr, "# saf"))
				formatType = CSequence::ePHDSAF;
			else
				formatType = CSequence::ePHD;
		}
			
		else if (strstr(searchStr, "MSF:"))
			formatType = CSequence::eMSF;

		// must be SAF or something

		// restore string
		if (endSearch)
			*endSearch = originalChar;
	}
	return formatType;

}


//_______________________________________________________________________________________
//
//	ParseAlignmentFromBuffer
//
//		Parse a text buffer and obtain the sequences
//
//_______________________________________________________________________________________
bool ParseAlignmentFromBuffer(char* buf, CParseMessageCallback& callback,
								CSequence& sequence, const CParseOptions& options)
{
	std::auto_ptr<CParse> parser;
	switch (sequence.FormatType)
	{
	case CSequence::eFASTA:
		{
			// ugly syntax because of lack of std compliance of auto_ptr in VC6 <memory>
			std::auto_ptr<CParse> temp(new CParseFASTA(callback));
			parser = temp;
		}
		break;
	case CSequence::ePHDSAF:
	case CSequence::eCLUSTAL:
		{
			// ugly syntax because of lack of std compliance of auto_ptr in VC6 <memory>
			std::auto_ptr<CParse> temp(new CParseCLUSTAL(callback));
			parser = temp;
		}
		break;
	case CSequence::ePHD:
	case CSequence::eDSC:
	case CSequence::eMSF:
		{
			// ugly syntax because of lack of std compliance of auto_ptr in VC6 <memory>
			std::auto_ptr<CParse> temp(new CParseMSF(callback,
														options.m_Ignore_COIL,
														options.m_IgnoreConfidences));
			parser = temp;
		}
		break;
	default:
		return false;
	}
	if (!parser->DoParse(buf, sequence))
		return false;

	if (sequence.FormatType == CSequence::ePHD || sequence.FormatType == CSequence::ePHDSAF)
	{
		CParsePHD parsePHD(callback, parser->GetLineNum(),
										options.m_Ignore_COIL,
										options.m_MergePHDPredictions,
										options.m_IgnoreConfidences);
		if (!parsePHD.DoParse(parser->GetDataStream(), sequence))
			return false;
	}
	return true;
}

//_______________________________________________________________________________________
//
//	CreateMultiSequences
//
//		Splits and copies the original alignment into multiple sets of sequences
//
//		Input:
//				indices		: indices into the original alignment
//								-1 denotes new sets
//				originalSeq	: Original alignment
//
//		Output:
//				vecSeq		: Sets of alignments
//
//		Logic:
//				Turn the indices into sets of indices into the original alignment
//				Then create the mulitsets by looking up the actual data
//
//_______________________________________________________________________________________
namespace{
	bool Contains(const std::vector<int>& vecIndices, int i)
	{
		return std::find(vecIndices.begin(), vecIndices.end(), i) != vecIndices.end();
	}
}
void CreateMultiSequences(std::vector<int>& indices, const CSequence originalSeq,
													std::vector<CSequence>& vecSeq)
{
	// indices [0] = beginning of first set should == -1
	assert(indices[0] == -1);
	if (indices.size() < 2)
		return;

	unsigned curSet = 0;
	std::vector<std::vector<int> > vecvecIndex;
	vecvecIndex.push_back(std::vector<int>());

	// get vector of indices corresponding to the sets of sequences
	for (unsigned i = 1; i < indices.size(); ++i)
	{
		int index = indices[i];
		// -1 == start new set of indices
		if (index < 0)
		{
			assert(index == -1);
			// only start new set if present set not empty
			if (vecvecIndex[curSet].size() > 0)
			{
				++curSet;
				vecvecIndex.push_back(std::vector<int>());
			}
			continue;
		}

		// ignore duplicates within a set . A sequence can only exist once per alignment
		if (index < (int)originalSeq.Count() &&
			!Contains(vecvecIndex[curSet], index))
			vecvecIndex[curSet].push_back(index);
	}


	// now turn indices into sets of sequences
	vecSeq.clear();
	for (unsigned i = 0; i < vecvecIndex.size(); ++i)
	{
		// break on (last) empty set
		if (!vecvecIndex[i].size())
		{
			// only empty sets should be at the end
			assert(i == vecvecIndex.size() - 1);
			break;
		}

		// create new sequence set
		vecSeq.push_back(CSequence());

		// lookup from indices
		for (unsigned j = 0; j < vecvecIndex[i].size(); ++j)
		{
			unsigned index = vecvecIndex[i][j];
			vecSeq[i].Data.push_back(originalSeq.Data[index]);
			vecSeq[i].Names.push_back(originalSeq.Names[index]);
		}
		vecSeq[i].Len = originalSeq.Len;
		assert(vecSeq[i].IS_VALID());
	}


}

//_______________________________________________________________________________________
//
//	OutputAnnotatedSequence
//
//		Description:
//			Takes the lists of amino acids and drives the OutputEngine to
//				produce annotated sequences.
//			Change the output type by passing in a different OutputEngine
//			Two overloaded versions are provided, the 2nd of which is optimized
//				for when the number of protein sequences < 32.
//
//		Inputs:
//			sequences			: Sets of aligned amino acid sequences
//			formatSequences		: To hold the format (indices) of corresponding sequences
//									at corresponding positions.
//			wrap				: Long alignments wrap around every "wrap" residues
//									To avoid wrapping, set wrap = sequence.Len
//			columnW				: Divide annotation into columns of columnW width
//									To avoid columns, set columnW = sequence.Len
//			OutputEngine		: Use this object to produce annotated output in
//									different file formats
//
//		Last modified: lg 11/04/2000
//
//_______________________________________________________________________________________
void OutputAnnotatedSequence(	const std::vector<CSequence>& sequences,
							    const std::vector<CSequence>& formatSequences,
							    unsigned wrap,
							    unsigned columnW,
							    CFormattedOutput& OutputEngine)
{
    // PRECONDITIONS
#ifndef NDEBUG
	assert(sequences.size());
	assert(sequences.size() == formatSequences.size());
	for (unsigned k = 0; k < sequences.size(); ++k)
		assert(sequences[k].Count() == formatSequences[k].Count());
	for (unsigned k = 0; k < sequences.size(); ++k)
		assert(sequences[k].IS_VALID());
	for (unsigned k = 0; k < sequences.size(); ++k)
		assert(formatSequences[k].IS_VALID());
	assert(wrap > 0 && wrap <= sequences[0].Len);
	assert(columnW > 0 && wrap <= sequences[0].Len);
#endif
    // PRECONDITIONS END


	for (unsigned beg = 0; beg < sequences[0].Len; beg += wrap)
	{
		unsigned end = std::min(sequences[0].Len, beg + wrap);
		for (unsigned k = 0; k < sequences.size(); ++k)
		{
			// add one and a half line sets
			if (k)
				OutputEngine.AddNewLineAndHalf();

			for (unsigned i = 0; i < sequences[k].Count(); ++i)
			{
				if (i)
					OutputEngine.AddNewLine();
				int currentFormat = 0;
				OutputEngine.SetFormat(0);
				OutputEngine.AddText(sequences[k].Names[i]);
				OutputEngine.AddText(' ');
				for (unsigned j = beg; j < end; ++j)
				{
					// space out columns
					if (!((j - beg) % columnW))
					{
						if (0 != currentFormat)
							OutputEngine.SetFormat(currentFormat = 0);
						OutputEngine.AddText(' ');
					}
					if (formatSequences[k].Data[i][j] != currentFormat)
	                    OutputEngine.SetFormat(currentFormat = formatSequences[k].Data[i][j]);

					OutputEngine.AddText(sequences[k].Data[i][j]);
				}
			}
			OutputEngine.SetFormat(0);

			// add paragraph to separate data from extra lines
			if (sequences[k].Count())
				OutputEngine.AddNewLine();

			// add extra lines, e.g. for consensus, 2ndary structure
			for (unsigned i = 0; i < sequences[k].ExtraCount(); ++i)
			{
				if (i)
					OutputEngine.AddNewLine();
				OutputEngine.AddText(sequences[k].ExtraNames[i]);
				OutputEngine.AddText(' ');
				for (unsigned j = beg; j < end; ++j)
				{
					// space out columns
					if (!((j - beg) % columnW))
						OutputEngine.AddText(' ');
					OutputEngine.AddText(sequences[k].ExtraData[i][j]);
				}
			}
		}

		// add one line between wraps
		OutputEngine.AddNewLine();
		OutputEngine.AddDoubleNewLines();
	}


}


